PNa <- seq(1,50,2)
R <- 8.314
Fa <- 96500
Temp <- 310
ghk <- function(PNa, PK, PCl, Na0, K0, Cl0, Nai, Ki, Cli){
Erev <- ((R*Temp)/Fa)*log((((PNa*Na0)+(PK*K0)+(PCl*Cli))/((PNa*Nai)+(PK*Ki)+(PCl*Cl0))))
return(Erev)
}
ghk(PNa,20,5,460,10,540,50,400,100)
## [1] -0.0594653967 -0.0441165074 -0.0345798147 -0.0276787293 -0.0222886434
## [6] -0.0178792788 -0.0141580032 -0.0109461228 -0.0081264758 -0.0056180724
## [11] -0.0033625875 -0.0013166234 0.0005529843 0.0022721489 0.0038615024
## [16] 0.0053377389 0.0067145554 0.0080033288 0.0092136107 0.0103534978
## [21] 0.0114299125 0.0124488190 0.0134153921 0.0143341503 0.0152090616
mempot <- ghk(PNa,20,5,460,10,540,50,400,100)
plot(PNa, mempot, xlab="Na+ permeability", ylab="membrane potential(V)", main="Membrane potential vs Na+ permeability",type="o")

P_Na <- seq(1,50,2)
mempot <- ghk(P_Na,20,5,460,10,540,50,400,100)
CURstate <- 0
endtime <- 100
state<-rep(0, endtime)
for(t in mempot){
if(t> -0.02){
p <- 0.8
q <- 0.8
} else {
p <- 0.65
q <- 0.03
}
for(t in 1:endtime){
randomNum<-runif(1) #runif(1) randomly generates a number between 0 and 1
if(CURstate==1){
if(randomNum < p){
CURstate <- 1
}else{
CURstate <- 0
}
}else{
if(randomNum< q){
CURstate<-1
}else{
CURstate<- 0
}
}
state[t] <- CURstate
}
mean(state)
plot(1:endtime, state, xlab="Time", ylab="Ion Channel state", main="ion channel state over time", type="o")
abline(h=mean(state))
legend("topleft",bty="n",legend=paste("mean",mean(state)), text.col=c(2))
}

























state
## [1] 0 1 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 1 0 0 1 1 1 0
## [38] 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 0 1
## [75] 0 1 1 0 1 1 1 0 1 1 1 1 1 1 0 1 1 1 0 1 1 1 1 0 1 1
#22. Overall, as the membrane potential gets more positive than the threshold value, the channel will be opened more often than it's closed, as when they are less than the threshold value -0.02 V, they will be closed more than they are open. For all the graphs, they retain one state more often than the other.